function upwind

%  Solves the advection equation using upwind
%       diff(u,t)+a*diff(u,x)=0   for xL < x < xr, 0 < t < tmax
%  where
%     u(x,0) = g(x)  and u = 0  at x=xL,xR  

% clear all previous variables and plots
clear *
clf

% set parameters
N=100;
M=72;
tmax=7;
xL=0;
xR=10;
a=1;

% pick time points (by picking the index)
itotal=4;
it(1)=1;
it(2)=round((M+1)/3);
it(3)=round(2*(M+1)/3);
it(4)=M+1;

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

% generate the points along the t-axis, t(1)=0 and t(M+1)=tmax
t=linspace(0,tmax,M+1);
k=t(2)-t(1);
lamda=k/h;

% calculate solutions
ue=exact(x,t,N+2,M+1,a);
u=fb(x,t,N+2,M+1,lamda);

% plot results
% get(gcf)
%set(gcf,'Position', [1290 314 560 725]);
plotsize(560,725)
for itt=1:itotal

	% plot solutions
	subplot(4,1,5-itt)
	plot(x,u(:,it(itt)),'-r')
	hold on
	plot(x,ue(:,it(itt)),'--k')
	
	% put in a directional arrow
	text(1.1+a*t(it(itt)),0.8,'\rightarrow','FontSize',18,'FontWeight','bold');
	
	% define axes and title used in plot
	ylabel('Solution','FontSize',14,'FontWeight','bold')

	% have MATLAB use certain plot options (all are optional)
	set(gca,'FontSize',14); 
	axis([xL xR -0.1 1.1]);
	box on
	
	% put in time value, label x-axis, and put in title
	if itt==1
		yt=0.95;
		xlabel('x-axis','FontSize',14,'FontWeight','bold')
	elseif itt==2
		yt=0.95;
	elseif itt==3
		yt=0.95;
	else
		yt=0.95;
		say2=['Upwind vs Exact Solution:  \lambda = ',num2str(lamda,'%3.3f'),'  (M = ',num2str(M),', N = ',num2str(N),')'];
		title(say2,'FontSize',14,'FontWeight','bold')
	end
	say=['t = ', num2str(t(it(itt)),'%3.2f')]; 
	text(8.5,yt,say,'FontSize',14,'FontWeight','bold')
	hold off
end;


% exact solution
function ue=exact(x,t,N,M,a)
for i=1:N
	ue(i,1)=g(x(i));
end;
for j=2:M
	for i=1:N
	  z=x(i)-a*t(j);
	  ue(i,j)=g(z);
	end;
end;

% F/t F/x method
function U=ff(x,t,N,M,lamda)
for i=1:N
	U(i,1)=g(x(i));
end;
for j=2:M
	U(1,j)=0;
	U(N,j)=0;
	for i=2:N-1
	  U(i,j)=U(i,j-1)-lamda*(U(i+1,j-1)-U(i,j-1));
	end;
end;

% F/t B/x method
function U=fb(x,t,N,M,lamda)
for i=1:N
	U(i,1)=g(x(i));
end;
for j=2:M
	U(1,j)=0;
	U(N,j)=0;
	for i=2:N-1
	  U(i,j)=U(i,j-1)-lamda*(U(i,j-1)-U(i-1,j-1));
	end;
end;

% subfunction g(x)
function q=g(x)
q=0;
if (x<1)&(x>0)
	q=1;
end;

% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);